package scu.maqiang.numeric;

import static java.lang.Math.*;

import java.util.Arrays;
import java.util.stream.IntStream;


/**
 * 
 * @author Administrator
 *
 */
public class IterSSolver extends Solver{
	
	public IterSSolver() {
	}
	
	public IterSSolver(SRMatrix mat) {
		int n = mat.getN();
		double[] diag = new double[n];
		CSRMatrix matU = new CSRMatrix(n);
		CSRMatrix matL = new CSRMatrix(n);
		CSRMatrix.extractDLU(mat, diag, matU, matL);
	}
	
	public IterSSolver(double eps, int maxIteration) {
		this.eps = eps;
		this.maxIteration = maxIteration;
	}

	/**
	 * �����ݶȷ����Գ��������Է�����(�Ѵ�����)
	 * @param A ϡ�����
	 * @param b �Ҷ���
	 * @param x0  ������ʼֵ, ͬʱ���������
	 * @return ��������, ��������������������������-1
	 */
	public void CG(SRMatrix A, double[] b, double[] x, int verbosity) {
		int n = A.getN();
		long begin = System.nanoTime();
		double[] Ap = A.mul(x);
		double[] r = MVO.add(b, -1.0, Ap);
		double[] p = (double[]) r.clone();
		double r2 = MVO.dot_product(r, r);
		int i;
		for (i = 0; i < maxIteration; i++) {
			if (sqrt(r2) < eps) {
				break;
			}
			//A.mul(p, Ap);
			//Ap = A.mul(p);
			//double temp1 = MVO.dot_product(Ap, p);
			double temp1 = 0.0;
			for(int k = 0; k < n; k++) {
				Ap[k] = 0.0;
				for(SREntry ele : A.mat[k]) {
					Ap[k] += ele.val * p[ele.col];
				}
				temp1 += Ap[k] * p[k];
			}
			double alpha = r2 / temp1;
			temp1 = 0;
			for (int k = 0; k < n; k++) {
				x[k] += alpha * p[k];
				r[k] -= alpha * Ap[k];
				temp1 += r[k] * r[k];
			}
			double beta = temp1 / r2;
			r2 = temp1;
			for (int k = 0; k < n; k++) {
				p[k] = r[k] + beta * p[k];
			}
			if (verbosity >= 2) {
				System.out.println();
				System.out.println("Iteration: " + (i + 1));
				System.out.printf("min: %20.10f\t max: %20.10f\n", MVO.min(x), MVO.max(x));
				System.out.printf("Error: %20.10f\n", sqrt(r2));
				System.out.println();
			}
		}
		if (verbosity >= 1) {
			System.out.println();
			System.out.println("Total iteration: " + i);
			System.out.printf("min: %20.10f\t max: %20.10f\n", MVO.min(x), MVO.max(x));
			long end = System.nanoTime();
			System.out.println("Solve time: " + (end - begin)/1.0e9 + " s");
			System.out.println();
		}
	}

	
	public void CG(SRMatrix A, double[] b, double[] x) {
		CG(A, b, x, 1);
	}
	
	public void CG(CSRMatrix A, double[] b, double[] x, int verbosity) {
		int n = A.getN();
		long begin = System.nanoTime();
		double[] Ap = A.mul(x);
		double[] r = MVO.add(b, -1.0, Ap);
		double[] p = (double[]) r.clone();
		double r2 = MVO.dot_product(r, r);
		double temp1, alpha, beta;
		int i, k;
		for (i = 0; i < maxIteration; i++) {
			if (sqrt(r2) < eps) {
				break;
			}
			A.mul(p, Ap);
			temp1 = MVO.dot_product(Ap, p);
			alpha = r2 / temp1;
			temp1 = 0;
			for (k = 0; k < n; k++) {
				x[k] += alpha * p[k];
				r[k] -= alpha * Ap[k];
				temp1 += r[k] * r[k];
			}
			beta = temp1 / r2;
			r2 = temp1;
			for (k = 0; k < n; k++) {
				p[k] = r[k] + beta * p[k];
			}
			if (verbosity >= 2) {
				System.out.println();
				System.out.println("Iteration: " + (i + 1));
				System.out.printf("min: %20.10f\t max: %20.10f\n", MVO.min(x), MVO.max(x));
				System.out.printf("Error: %20.10f\n", sqrt(r2));
				System.out.println();
			}
		}
		if (verbosity >= 1) {
			System.out.println();
			System.out.println("Total iteration: " + i);
			System.out.printf("min: %20.10f\t max: %20.10f\n", MVO.min(x), MVO.max(x));
			long end = System.nanoTime();
			System.out.println("Solve time: " + (end - begin)/1.0e9 + " s");
			System.out.println();
		}
	}
	
	/**
	 * Jacobi������������Է�����
	 * @param A ϡ�����
	 * @param b �Ҷ���
	 * @param x0 ������ʼֵ, ͬʱ���������
	 * @return ��������, ��������������������������-1
	 */
	public int Jacobi(SRMatrix A, double[] b, double[] x0) {
		int n = A.getN();
		double[] x1 = x0.clone();
		double err = eps + 1.0;
		double a00 = 0.0;
		int i;
		for (i = 0; i < maxIteration; i++) {
			if(err < Constants.Er) {
				break;
			}
			for(int j = 0; j < n; j++) {
				double temp = 0.0;
				for (SREntry ele: A.mat[j]) {
					if(ele.col != j) {
						temp += ele.val * x0[ele.col];
					}else {
						a00 = ele.val;
					}
				}
				x1[j] = (b[j] - temp)/a00;
			}
			err = MVO.distance(x0, x1);
			System.arraycopy(x1, 0, x0, 0, n);
		}
		if(i == maxIteration) {
			System.out.println("Iteration Has Reached Maximum Time! ");
			return -1;
		}
		System.out.println("Iteration Time " + i);
		return i;
	}
	
	/**
	 * GaussSeidel������������Է�����
	 * @param A ϡ�����
	 * @param b �Ҷ���
	 * @param x0 ������ʼֵ, ͬʱ���������
	 * @return ��������, ��������������������������-1
	 */
	public int GaussSeidel(SRMatrix A, double[] b, double[] x0) {
		A.sort();
		int n = A.getN();
		double err = eps + 1.0;
		int i;
		for(i = 0; i < maxIteration; i++) {
			if(err < Constants.Er) {
				break;
			}
			double L2Norm = MVO.L2Norm(x0);
			for(int j = 0; j < n; j++) {
				double temp = 0.0;
				double a00 = 0.0;
				for(SREntry ele: A.mat[j]) {
					if(ele.col != j) {
						temp += x0[ele.col] * ele.val;
					}else {
						a00 = ele.val;
					}
				}
				x0[j] = (b[j] - temp)/a00;
			}
			err = Math.abs((MVO.L2Norm(x0) - L2Norm));
		}
		if(i == maxIteration) {
			System.out.println("Iteration Has Reached Maximum Time! ");
			return -1;
		}
		System.out.println("Iteration Time " + i);
		return i;
	}
	
	/**
	 * ��γ��ɳڵ�����(Successive Over Relaxation(SOR))������Է�����
	 * @param A ϡ�����
	 * @param b �Ҷ���
	 * @param x0 ������ʼֵ, ͬʱ���������
	 * @param omega �ɳ�����, 0 < w < 2
	 * @return ��������, ��������������������������-1
	 */
	public int SOR(SRMatrix A, double[] b, double[] x0, double omega) {
		A.sort();
		int n = A.getN();
		double err = eps + 1.0;
		int i;
		for(i = 0; i < maxIteration; i++) {
			if(err < Constants.Er) {
				break;
			}
			double L2Norm = MVO.L2Norm(x0);
			for(int j = 0; j < n; j++) {
				double temp = 0.0;
				double a00 = 0.0;
				for(SREntry ele: A.mat[j]) {
					if(ele.col != j) {
						temp += x0[ele.col] * ele.val;
					}else {
						a00 = ele.val;
					}
				}
				x0[j] = (1-omega)*x0[j] + omega*(b[j] - temp)/a00;
			}
			err = Math.abs((MVO.L2Norm(x0) - L2Norm));
		}
		if(i == maxIteration) {
			System.out.println("Iteration Has Reached Maximum Time! ");
			return -1;
		}
		System.out.println("Iteration Time " + i);
		return i;
	}
	
	/**
	 * �Գ���γ��ɳڵ�����(Symmetric Successive Over Relaxation(SSOR))������Է�����
	 * @param A ϡ�����
	 * @param b �Ҷ���
	 * @param x0 ������ʼֵ, ͬʱ���������
	 * @param omega �ɳ�����, 0 < w < 2
	 * @return ��������, ��������������������������-1
	 */
	public int SSOR(SRMatrix A, double[] b, double[] x0, double omega) {
		A.sort();
		int n = A.getN();
		double err = eps + 1.0;
		int i;
		for(i = 0; i < maxIteration; i++) {
			if(err < Constants.Er) {
				break;
			}
			double L2Norm = MVO.L2Norm(x0);
			for(int j = 0; j < n; j++) {
				double temp = 0.0;
				double a00 = 0.0;
				for(SREntry ele: A.mat[j]) {
					if(ele.col != j) {
						temp += x0[ele.col] * ele.val;
					}else {
						a00 = ele.val;
					}
				}
				x0[j] = (1-omega)*x0[j] + omega*(b[j] - temp)/a00;
			}
			
			for(int j = n - 1; j >= 0; j--) {
				double temp = 0.0;
				double a00 = 0.0;
				for(SREntry ele: A.mat[j]) {
					if(ele.col != j) {
						temp += x0[ele.col] * ele.val;
					}else {
						a00 = ele.val;
					}
				} 
				x0[j] = (1-omega)*x0[j] + omega*(b[j] - temp)/a00;
			}
			err = Math.abs((MVO.L2Norm(x0) - L2Norm));
		}
		if(i == maxIteration) {
			System.out.println("Iteration Has Reached Maximum Time! ");
			return -1;
		}
		System.out.println("Iteration Time " + i);
		return i;
	}
	
	/**
	 * Ԥ�������ݶȷ����Գ��������Է�����, ����������ĶԽ�Ԫ�ؽ���Ԥ����(�Ѵ�����)
	 * @param A ϡ�����
	 * @param b �Ҷ���
	 * @param x0  ������ʼֵ, ͬʱ���������
	 * @return ��������, ��������������������������-1
	 */
	public int PCGM(SRMatrix A, double[] b, double[] x0) {
		int n = A.getN();
		
		double[] r0 = MVO.add(b, -1.0, A.mul(x0));
		
		double[] z0 = new double[n];
		double[] M = A.diag();
		for (int i = 0; i < z0.length; i++) {
			z0[i] = r0[i] / M[i];
		}
		
		double[] z1 = z0.clone();
		double[] p0 = z0.clone();
		
		int i;
		for (i = 0; i < maxIteration; i++) {
			if (MVO.L2Norm(r0) < eps) {
				break;
			}
			double temp1 = 0.0;
			double temp2 = 0.0;
			for (int j = 0; j < n; j++) {
				z1[j] = 0.0;
				for (SREntry ele : A.mat[j]) {
					z1[j] += p0[ele.col] * ele.val;
				}
				temp1 += z0[j] * r0[j];
				temp2 += z1[j] * p0[j];
			}
			
			double alpha = temp1 / temp2; 
			temp2 = temp1;
			temp1 = 0.0;
			for (int j = 0; j < n; j++) {
				x0[j] += alpha * p0[j];
				r0[j] -= alpha * z1[j];
				z1[j] = r0[j] / M[j];
				temp1 += z1[j] * r0[j];
			}
			double beta = temp1 / temp2; 
			for (int j = 0; j < n; j++) {
				p0[j] = z1[j] + beta * p0[j];
				z0[j] = z1[j];
			}
		}
		if(i == maxIteration) {
			System.out.println("Iteration Has Reached Maximum Time! ");
			return -1;
		}
		System.out.println("Iteration Time " + i);
		return i;
	}
	
	/**
	 * �Գ���γ��ɳ�Ԥ�������ݶȷ����Գ��������Է�����(�Ѵ�����)
	 * @param A �Գ�����ϡ�����
	 * @param b �Ҷ���
	 * @param x0  ������ʼֵ, ͬʱ���������
	 * @param omega �ɳ�����
	 * @return ��������, ��������������������������-1
	 */
	public int PCGSSOR(SRMatrix A, double[] b, double[] x0, double omega, int verbosity) {
		long begin = System.nanoTime();
		int n = A.getN();
		double[] r0 = MVO.add(b, -1.0, A.mul(x0));
		
		double[] z0 = new double[n];
		double[] diag = A.diag();
		SSORTreat(A, diag, r0, z0, omega);
		//System.out.println(Arrays.toString(z0));
		double[] z1 = z0.clone();
		double[] p0 = z0.clone();
		
		int i, j;
		for (i = 0; i < maxIteration; i++) {
			if (MVO.L2Norm(r0) < eps) {
				break;
			}
			double temp1 = 0.0;
			double temp2 = 0.0;
			for (j = 0; j < n; j++) {
				z1[j] = 0.0;
				for (SREntry ele : A.mat[j]) {
					z1[j] += p0[ele.col] * ele.val;
				}
				temp1 += z0[j] * r0[j];
				temp2 += z1[j] * p0[j];
			}
			
			//System.out.println(temp1 +"\t" + temp2);
			double alpha = temp1 / temp2; 
			for (j = 0; j < n; j++) {
				x0[j] += alpha * p0[j];
				r0[j] -= alpha * z1[j];
			}
			
			//System.out.println("r0 = " + MVO.toString(r0));
			SSORTreat(A, diag, r0, z1, omega);
			
			temp2 = temp1;
			temp1 = 0.0;
			for (j = 0; j < n; j++) {
				temp1 += z1[j] * r0[j];
			}
			double beta = temp1 / temp2; 
			for (j = 0; j < n; j++) {
				p0[j] = z1[j] + beta * p0[j];
				z0[j] = z1[j];
			}
			if (verbosity >= 2) {
				System.out.println();
				System.out.println("Iteration: " + (i + 1));
				System.out.printf("min: %20.10f\t max: %20.10f\n", MVO.min(x0), MVO.max(x0));
				System.out.printf("Error: %20.10f\n", MVO.L2Norm(r0));
				System.out.println();
			}
			//System.arraycopy(r1, 0, r0, 0, n);
			//System.arraycopy(z1, 0, z0, 0, n);
		}
		
		if (verbosity >= 1) {
			System.out.println();
			System.out.println("Total iteration: " + i);
			System.out.printf("min: %20.10f\t max: %20.10f\n", MVO.min(x0), MVO.max(x0));
			long end = System.nanoTime();
			System.out.println("Solve time: " + (end - begin)/1.0e9 + " s");
			System.out.println();
		}
		if(i == maxIteration) {
			System.out.println("Iteration Has Reached Maximum Time! ");
			return -1;
		}
		return i;
	}

	public int PCGSSOR_2(SRMatrix A, double[] b, double[] x0, double omega, int verbosity) {
		int n = A.getN();
		double[] diag = new double[n];
		CSRMatrix matU = new CSRMatrix(n);
		CSRMatrix matL = new CSRMatrix(n);
		CSRMatrix.extractDLU(A, diag, matU, matL);
		
		long begin = System.nanoTime();
		double[] r0 = MVO.add(b, -1.0, matU.mul(x0));
		double[] z0 = new double[n];
		CSRMatrix.SSORTreat(diag, matU, matL, r0, z0, omega);
		
		double[] z1 = z0.clone();
		double[] p0 = z0.clone();
		
		int i, j;
		for (i = 0; i < maxIteration; i++) {
			if (MVO.L2Norm(r0) < eps) {
				break;
			}
			double temp1 = 0.0;
			double temp2 = 0.0;
			
			matU.mul(p0, z1);
			for (j = 0; j < n; j++) {
				temp1 += z0[j] * r0[j];
				temp2 += z1[j] * p0[j];
			}
			double alpha = temp1 / temp2; 
			for (j = 0; j < n; j++) {
				x0[j] += alpha * p0[j];
				r0[j] -= alpha * z1[j];
			}
			
			//System.out.println("r0 = " + MVO.toString(r0));
			//SSORTreat(A, diag, r0, z1, omega);
			CSRMatrix.SSORTreat(diag, matU, matL, r0, z1, omega);
			
			temp2 = temp1;
			temp1 = 0.0;
			for (j = 0; j < n; j++) {
				temp1 += z1[j] * r0[j];
			}
			double beta = temp1 / temp2; 
			for (j = 0; j < n; j++) {
				p0[j] = z1[j] + beta * p0[j];
				z0[j] = z1[j];
			}
			if (verbosity >= 2) {
				System.out.println();
				System.out.println("Iteration: " + (i + 1));
				System.out.printf("min: %20.10f\t max: %20.10f\n", MVO.min(x0), MVO.max(x0));
				System.out.printf("Error: %20.10f\n", MVO.L2Norm(r0));
				System.out.println();
			}
			//System.arraycopy(r1, 0, r0, 0, n);
			//System.arraycopy(z1, 0, z0, 0, n);
		}
		
		if (verbosity >= 1) {
			System.out.println();
			System.out.println("Total iteration: " + i);
			System.out.printf("min: %20.10f\t max: %20.10f\n", MVO.min(x0), MVO.max(x0));
			long end = System.nanoTime();
			System.out.println("Solve time: " + (end - begin)/1.0e9 + " s");
			System.out.println();
		}
		if(i == maxIteration) {
			System.out.println("Iteration Has Reached Maximum Time! ");
			return -1;
		}
		System.out.println("Iteration Time " + i);
		return i;
	}
	/**
	 * ���öԳ���γ��ɳڵ�������Ԥ�������
	 * @param A �Գ�����ϡ�����
	 * @param diag ����Խ���Ԫ��
	 * @param r0 ����
	 * @param z0 
	 * @param omega �ɳ�����
	 */
	private void SSORTreat(SRMatrix A, double[] diag, double[] r0, double[] z0, double omega) {
		int n = A.getN();
		
		z0[0] = r0[0] / diag[0];
		for (int i = 1; i < n; i++) {
			double temp = 0.0;
			for (SREntry ele : A.mat[i]) {
				if(ele.col < i) {
					temp += omega * ele.val * z0[ele.col];
				}
			}
			//double ii = i;
			//double temp = A.mat[i].stream().filter(ele -> ele.col < ii).mapToDouble(ele -> omega * ele.val * z0[ele.col]).sum();
			z0[i] = (r0[i] - temp) / diag[i];
		}
		//System.out.println("z0 = " + MVO.toString(z0));
		for(int i = 0; i < n; i++) {
			z0[i] = z0[i] * diag[i];
		}
		z0[n - 1] = z0[n - 1] / diag[n - 1];
		for (int i = n - 2; i >= 0; i--) {
			double temp = 0.0;
			for(SREntry ele : A.mat[i]) {
				if(ele.col > i) {
					temp += omega * ele.val * z0[ele.col];
				}
			}
			z0[i] = (z0[i] - temp) / diag[i];
		}
	}
	
	/**
	 * �����ݶȷ��������(CGNE)�����һ�����Է�����(�Ѵ�����)
	 * @param A һ��ϡ�����
	 * @param b �Ҷ���
	 * @param x0 ������ʼֵ�������ص�����
	 * @return ��������, ��������������������������-1
	 */
	public int CGNE(SRMatrix A, double[] b, double[] x0) {
		int n = b.length;
		double[] r0 = MVO.add(b, -1.0, A.mul(x0));
		double[] p0 = A.mulT(r0);
		double[] Ap0 = new double[n];
		int i;
		for (i = 0; i < maxIteration; i++) {
			if (MVO.L2Norm(r0) < eps) {
				break;
			}
			double temp1 = 0.0;
			double temp2 = 0.0;
			for (int j = 0; j < n; j++) {
				Ap0[j] = 0.0;
				for (SREntry ele : A.mat[j]) {
					Ap0[j] += p0[ele.col] * ele.val;
				}
				temp1 += r0[j] * r0[j];
				temp2 += p0[j] * p0[j];
			}
			double alpha = temp1 / temp2 ; 
			temp2 = temp1;
			temp1 = 0.0;
			for (int j = 0; j < n; j++) {
				x0[j] += alpha * p0[j];
				r0[j] -= alpha * Ap0[j];
				temp1 += r0[j] * r0[j];
			}
			double beta = temp1 / temp2;
			Arrays.fill(Ap0, 0.0);
			for (int j = 0; j < n; j++) {
				for (SREntry ele : A.mat[j]) {
					Ap0[ele.col] += r0[j] * ele.val;
				}
				p0[j] = Ap0[j] + beta * p0[j];
			}
			System.out.println("Iteration " + (i + 1) + ": min\t"
					+ MVO.min(x0) + "\tmax\t" + MVO.max(x0));
		}
		if (i == maxIteration) {
			System.out.println("Iteration Has Reached Maximum Time");
			return -1;
		}
		System.out.println("Iteration Time " + i);
		return i;
	}
	
	/**
	 * �����ݶȷ����̲���(CGNR)�����һ�����Է�����(�Ѵ�����)
	 * @param A һ�������ϡ�����
	 * @param b �Ҷ���
	 * @param x0 ������ʼֵ�������ص�����
	 * @return ��������, ��������������������������-1
	 */
	public int CGNR(SRMatrix A, double[] b, double[] x0) {
		int n = b.length;
		double[] r0 = MVO.add(b, -1.0, A.mul(x0));
		double[] p0 = A.mulT(r0);
		double[] Ap0 = new double[n];
		int i;
		for (i = 0; i < maxIteration; i++) {
			if (MVO.L2Norm(r0) < eps) {
				break;
			}
			double temp1 = 0.0;
			double temp2 = 0.0;
			Arrays.fill(Ap0, 0.0);
			for (int j = 0; j < n; j++) {
				for (SREntry ele : A.mat[j]) {
					Ap0[ele.col] += r0[j] * ele.val;
				}
				temp1 += Ap0[j] * Ap0[j];
				
				Ap0[j] = 0.0;
				for (SREntry ele : A.mat[j]) {
					Ap0[j] += p0[ele.col] * ele.val;
				}
				temp2 += Ap0[j] * Ap0[j];
			}
			
			double alpha = temp1 / temp2;
			for (int j = 0; j < n; j++) {
				x0[j] += alpha * p0[j];
				r0[j] -= alpha * Ap0[j];
			}
			
			temp2 = temp1;
			temp1 = 0.0;
			Arrays.fill(Ap0, 0.0);
			for (int j = 0; j < n; j++) {
				for (SREntry ele : A.mat[j]) {
					Ap0[ele.col] += r0[j] * ele.val;
				}
				temp1 += Ap0[j] * Ap0[j];
			}
			double beta = temp1 / temp2;
			for (int j = 0; j < n; j++) {
				p0[j] = Ap0[j] + beta * p0[j];
			}

			System.out.println("Iteration " + (i + 1) + ": min\t"
					+ MVO.min(x0) + "\tmax\t" + MVO.max(x0));
		}
		if (i == maxIteration) {
			System.out.println("Iteration Has Reached Maximum Time");
			return -1;
		}
		System.out.println("Iteration Time " + i);
		return i;
	}
	
	
	/**
	 * �����ݶ�ƽ����(CGS)���һ�����Է�����
	 * @param A һ�������ϡ�����
	 * @param b �Ҷ���
	 * @param x0 ������ʼֵ�������ص�����
	 * @return ��������, ��������������������������-1
	 */
	public int CGS(SRMatrix A, double[] b, double[] x0) {
		int n = A.getN();
		double[] r0 = MVO.add(b, -1.0, A.mul(x0));
		double[] u0 = r0.clone();
		double[] z0 = MVO.add(r0, 1.0);
		double[] p0 = r0.clone();
		double[] q0 = new double[n];
		double[] v0 = A.mul(p0);
		//double alpha = 0.0;
		//double beta = 0.0;
		int k;
		for (k = 0; k < maxIteration; k++) {
			if (MVO.L2Norm(r0) < eps) {
				break;
			}
			
			double temp1 = 0.0;
			double temp2 = 0.0;
			for(int j = 0; j < n; j++) {
				temp1 += r0[j] * z0[j];
				temp2 += v0[j] * z0[j];
			}
			double alpha = temp1 / temp2;
//			for(int i = 0; i < n; i++) {
//				q0[i] = u0[i] - alpha * v0[i];
//			}
//			double[] uq = MVO.add(u0, q0);
//			double[] Auq = A.mul(uq);
//			
//			temp2 = temp1;
//			temp1 = 0.0;
//			for(int i = 0; i < n; i++) {
//				x0[i] += alpha * uq[i];
//				r0[i] -= alpha * Auq[i];
//				temp1 += r0[i] * z0[i];
//			}
			for(int i = 0; i < n; i++) {
				q0[i] = u0[i] - alpha * v0[i];
				v0[i] = u0[i] + q0[i];

			}
			
//			System.out.println("v0 = " + MVO.toString(v0));
			//double[] uq = MVO.add(u0, q0);
			//double[] Auq = A.mul(v0);
			temp2 = temp1;
			temp1 = 0.0;
			for(int i = 0; i < n; i++) {
				u0[i] = 0.0;
				for (SREntry ele : A.mat[i]) {
					u0[i] += v0[ele.col] * ele.val;
				}
				x0[i] += alpha * v0[i];
				r0[i] -= alpha * u0[i];
				temp1 += r0[i] * z0[i];
			}
			//u0 = A.mul(v0);
			
			
//			for(int i = 0; i < n; i++) {
//				x0[i] += alpha * v0[i];
//				r0[i] -= alpha * u0[i];
//				temp1 += r0[i] * z0[i];
//			}
			
			double beta = temp1 / temp2;
			for(int i = 0; i < n; i++) {
				u0[i] = r0[i] + beta * q0[i];
				p0[i] = u0[i] + beta *(q0[i] + beta * p0[i]);

			}
			for(int i = 0; i < n; i++) {
				v0[i] = 0.0;
				for (SREntry ele : A.mat[i]) {
					v0[i] += p0[ele.col] * ele.val;
				}
			}
			//v0 = A.mul(p0);
//			System.out.println("Iteration " + (k + 1) + ": min\t"
//					+ MVO.min(x0) + "\tmax\t" + MVO.max(x0));
		}
		if(k == maxIteration) {
			System.out.println("Iteration Has Reached Maximum Time! ");
			return -1;
		}
		System.out.println("Iteration Time " + k);
		return k;
	}
	
	/**
	 * ˫�����ݶȷ�(BiCG)���һ�����Է�����
	 * @param A һ�������ϡ�����
	 * @param b �Ҷ���
	 * @param x0 ������ʼֵ�������ص�����
	 * @return ��������, ��������������������������-1
	 */
	public int BiCG(SRMatrix A, double[] b, double[] x0) {
		int n = A.getN();
		double[] r0 = MVO.add(b, -1.0, A.mul(x0));
		double[] r1 = r0.clone();
		double[] p0 = r0.clone();
		double[] z0 = r0.clone();
		double[] z1 = new double[n];
		double[] pt = z0.clone();
		int i;
		for(i = 0; i < maxIteration; i++) {
			if (MVO.L2Norm(r0) < eps) {
				break;
			}
			double[] Ap0 = A.mul(p0);
			double alpha = MVO.dot_product(r0, z0)/MVO.dot_product(Ap0, pt);
			double[] Atpt = A.mulT(pt);
			for(int j = 0; j < n; j++) {
				x0[j] = x0[j] + alpha * p0[j];
				r1[j] = r0[j] - alpha * Ap0[j];
				z1[j] = z0[j] - alpha * Atpt[j];
			}
			System.out.println("z1 = " + MVO.toString(z1));
			double beta = MVO.dot_product(r1, z1)/MVO.dot_product(r0, z0);
			System.out.println(beta);
			for(int j = 0; j < n; j++) {
				p0[j] = r1[j] + beta * p0[j];
				pt[j] = z1[j] + beta * pt[j];
				z0[j] = z1[j];
				r0[j] = r1[j];
			}
			System.out.println("Iteration " + (i + 1) + ": min\t"
					+ MVO.min(x0) + "\tmax\t" + MVO.max(x0));
		}
		if (i == maxIteration) {
			System.out.println("Iteration Has Reached Maximum Time");
			return -1;
		}
		System.out.println("Iteration Time " + i);
		return i;
	}
	
	/**
	 * ˫�����ݶ��ȶ�������������Է�����
	 * @param A һ�������ϡ�����
	 * @param b �Ҷ���
	 * @param x0 ������ʼֵ�������ص�����
	 * @return ��������, ��������������������������-1
	 */
	public int BiCGSTAB(SRMatrix A, double[] b, double[] x0) {
		int n = A.getN();
		double[] r0 = MVO.add(b, -1.0, A.mul(x0));
		double[] p0 = r0.clone();
		double[] z0 = r0.clone();
		int i;
		for(i = 0; i < maxIteration; i++) {
			if(MVO.L2Norm(r0) < Constants.Er) {
				break;
			}
			System.out.println("Error = " + MVO.L2Norm(r0));
			double[] Ap0 = A.mul(p0);
			double temp = MVO.dot_product(r0, z0);
			double alpha = temp/ MVO.dot_product(Ap0, z0);
			for(int j = 0; j < n; j++) {
				x0[j] = x0[j] + alpha * p0[j];
				r0[j] = r0[j] - alpha * Ap0[j];
			}
			double[] Ar0 = A.mul(r0);
			double omega = MVO.dot_product(r0, Ar0)/MVO.dot_product(Ar0, Ar0);
			for(int j = 0; j < n; j++) {
				x0[j] = x0[j] + omega * r0[j];
				r0[j] = r0[j] - omega * Ar0[j];
			}
			double beta = alpha /omega * MVO.dot_product(r0, z0) / temp;
			for(int j = 0; j < n; j++) {
				p0[j] = r0[j] + beta *(p0[j] - omega * Ap0[j]);
			}
			
//			System.out.println("Iteration " + (i + 1) + ": min\t"
//					+ MVO.min(x0) + "\tmax\t" + MVO.max(x0));
		}
		
		if (i == maxIteration) {
			System.out.println("Iteration Has Reached Maximum Time");
			return -1;
		}
		System.out.println("Iteration Time " + i);
		return i;
	}
	
	/**
	 * �Գ���γ��ɳ�Ԥ�������ݶ�ƽ����(CGS)���һ�����Է�����
	 * @param A һ�������ϡ�����
	 * @param b �Ҷ���
	 * @param x0 ������ʼֵ�������ص�����
	 * @return ��������, ��������������������������-1
	 */
	public int PCGSSSOR(SRMatrix A, double[] b, double[] x0, double omega) {
		int n = A.getN();
		double[] r0 = MVO.add(b, -1.0, A.mul(x0));
		
		double[] z0 = new double[n];
		double[] diag = A.diag();
		SSORTreat(A, diag, r0, z0, omega);
		
		double[] u0 = z0.clone();
		double[] p0 = z0.clone();
		double[] q0 = new double[n];
		double[] Ap0 = A.mul(p0);
		double[] v0 = new double[n];
		
		double[] rr0 = MVO.add(z0, 1.0);
		
		SSORTreat(A, diag, Ap0, v0, omega);
		//double alpha = 0.0;
		//double beta = 0.0;
		int k;
		for (k = 0; k < maxIteration; k++) {
			if (MVO.L2Norm(r0) < eps) {
				break;
			}
			System.out.println("Error = " + MVO.L2Norm(r0));
			double temp1 = 0.0;
			double temp2 = 0.0;
			for(int j = 0; j < n; j++) {
				temp1 += z0[j] * rr0[j];
				temp2 += v0[j] * rr0[j];
			}
			double alpha = temp1 / temp2;
//			for(int i = 0; i < n; i++) {
//				q0[i] = u0[i] - alpha * v0[i];
//			}
//			double[] uq = MVO.add(u0, q0);
//			double[] Auq = A.mul(uq);
//			
//			temp2 = temp1;
//			temp1 = 0.0;
//			for(int i = 0; i < n; i++) {
//				x0[i] += alpha * uq[i];
//				r0[i] -= alpha * Auq[i];
//				temp1 += r0[i] * z0[i];
//			}
			for(int i = 0; i < n; i++) {
				q0[i] = u0[i] - alpha * v0[i];
				v0[i] = u0[i] + q0[i];
			}
			//double[] uq = MVO.add(u0, q0);
			//double[] Auq = A.mul(v0);
			
			
			for(int i = 0; i < n; i++) {
				u0[i] = 0.0;
				for (SREntry ele : A.mat[i]) {
					u0[i] += v0[ele.col] * ele.val;
				}
				x0[i] += alpha * v0[i];
				r0[i] -= alpha * u0[i];
				//temp1 += r0[i] * rr0[i];
			}
			//u0 = A.mul(v0);
			
			SSORTreat(A, diag, r0, z0, omega);
			
			temp2 = temp1;
			temp1 = 0.0;
			for(int i = 0; i < n; i++) {
				temp1 += z0[i] * rr0[i];
			}
//			for(int i = 0; i < n; i++) {
//				x0[i] += alpha * v0[i];
//				r0[i] -= alpha * u0[i];
//				temp1 += r0[i] * z0[i];
//			}
			
			double beta = temp1 / temp2;
			for(int i = 0; i < n; i++) {
				u0[i] = z0[i] + beta * q0[i];
				p0[i] = u0[i] + beta *(q0[i] + beta * p0[i]);

			}
			for(int i = 0; i < n; i++) {
				Ap0[i] = 0.0;
				for (SREntry ele : A.mat[i]) {
					Ap0[i] += p0[ele.col] * ele.val;
				}
			}
			//Ap0 = A.mul(p0);
			SSORTreat(A, diag, Ap0, v0, omega);
			//v0 = A.mul(p0);
//			System.out.println("Iteration " + (k + 1) + ": min\t"
//					+ MVO.min(x0) + "\tmax\t" + MVO.max(x0));
		}
		if(k == maxIteration) {
			System.out.println("Iteration Has Reached Maximum Time! ");
			return -1;
		}
		System.out.println("Iteration Time " + k);
		return k;
	}
	
	public void GMRES(SRMatrix A, double[] b, double[] x, int m) {
		int n = A.getM();
		ParamCheck.checkEqual(n, b.length);
		ParamCheck.checkEqual(n, x.length);
		double[] r = MVO.add(b, -1.0, A.mul(x));
		double[][] H = new double[m + 1][m];
		double[][] Q = new double[n][m + 1];
		double[] x0 = x.clone();
		double[] f = new double[m + 1];
		f[0] = MVO.L2Norm(r);
		for (int i = 0; i < n; i++) {
			Q[i][0] = r[i] / f[0];
		}
		for (int k = 0; k < m; k++) {
			double[] colK = MVO.col(Q, k);
			double[] y = A.mul(colK);
			for (int j = 0; j <= k; j++) {
				double[] colJ = MVO.col(Q, j);
				H[j][k] = MVO.dot_product(colJ, y);
				for (int t = 0; t < n; t++) {
					y[t] -= H[j][k] * colJ[t];
				}
			}
			H[k + 1][k] = MVO.L2Norm(y);
			for (int t = 0; t < n; t++) {
				Q[t][k + 1] = y[t] / H[k + 1][k];
			}
			double[] c = new double[k + 1];
			double error = IterDSolver.MLS_QR_HH(MVO.subMatrix(H, 0, k + 1, 0, k), MVO.subVector(f, 0, k + 1), c);
			MVO.add(x0, 1.0, MVO.matmul(MVO.subMatrix(Q, 0, n - 1, 0, k), c), x);
			System.out.println("Error = " + error);
			if (Math.abs(error) < Constants.Er) {
				break;
			}
		}

	}
	private CSRMatrix matU;
	private CSRMatrix matL;
	private double[] diag;
}

